Analysis of sinusoidal post-buckling deformation of horizontal coiled tubing with initial residual bending

In the process of horizontal well construction, the working environment in the well is very complicated. It is difficult to determine the actual deformation and force of the coiled tubing in the wellbore from the load change at the wellhead along. In addition, the coiled tubing wrapped around the reel and entering the wellbore through the injection head may cause an initial residual bending. In this paper, the mathematical model of the coiled tubing under axial load during buckling deformation is established, and the analytical solution with the time term is derived by reasonable simplification of the model. The effect of residual bending on the post-sinusoidal buckling of the coiled tubing in the horizontal well was investigated using the separation variable method. The effect of residual bending variations on the well wall contact forces is analyzed. The research shows that the initial parameter Г is the main factor influencing the post-sinusoidal buckling of the coiled tubing with residual bending. The larger the parameter value Г, the greater the effect of the residual bending on the post-sinusoidal buckling deformation of the coiled tubing, and the earlier the critical point of the mixed sinusoidal-helical buckling of the coiled tubing appears. The compression velocity of the coiled tubing has a significant effect on the well wall contact force. The faster the compression, the greater the contact force. By introducing the time term and the separation variable, this paper provides a new method and theoretical basis for further study of the process of entering sinusoidal-helical buckling of the coiled tubing with initial residual bending.


Introduction
Coiled tubing (CT) has the advantages of fast and efficient operation, low formation damage, and low cost, and is now widely used in well workovers, cleanups, sidetracks, and even drilling operations.However, due to the unique characteristics of the CT, buckling behavior usually occurs when working in complex downhole environments [1,2], causing significant increases in downhole frictional resistance and torque, even locking, or causing failure of the tubular string connections [3][4][5].Since Rubinsky pioneered study of the deformation of strings during operation and provided the differential equations for strings buckling for the first time, many scholars have also conducted research on the buckling of the tubular string [6].Mitchell derived the differential equations for helical buckling based on the principle of virtual work, and used the Lunger-Kutta method to find the analytical solution [7].Paslay and Dawson obtained the lateral stability conditions for inclined straight wells [8,9].Chen et al. developed a spiral buckling stability criterion [10], and Mitchell investigated the effect of wellbore curvature on the buckling stability of the tubular string [11].Huang and Pattillo, Hishida, and other scholars have conducted in-depth research on helical buckling [12,13].Sun et al. proposed an improved helical buckling model that includes both helical and non-helical segments to reduce the error [14].
With further research, many scholars have found that the deformation of the tubular string in complex wellbores is affected by various conditions [15][16][17].Mitchell, Graphics, and Miska, Gao, and Di et al. studied the deformation state of the tubular string with joints during operation [18,19].Liu et al. studied the sinusoidal buckling and helical buckling of the tubular string in the wellbore by using the inextensible and extensible rod theories [20], and pointed out that sinusoidal buckling deformation occurs firstly in the tubular string compression process, and then it enters the mixed sinusoidal-helical buckling deformation, and finally becomes the stabilized helical buckling deformation.However, so far there is no consistent mathematical model to represent the conversion from sinusoidal to helical buckling of the tubular string [21].Gao studied the velocity variation during deformation of the tubular string [22].Gao also pointed out that the velocity of movement of the tubular string has little effect on the critical load of the tubular string buckling [23,24].However, the important feature of the tubular string during the conversion morphology is the movement pattern of the micro segments of the tubular string.
In production, the coiled tubing is wound on a reel.The plastic deformation that occurs after a long period of time remains within the metal structure of the coiled tubing.Even if the CT is straightened, the plastic deformation remains for a long time, and it is called "residual bending".The CT is lowered into the wellbore by the rotation of the reel.Without the centralizer, the CT entering the wellbore may have an initial curvature.Qiu and Miska, and Zhu et al. have studied the buckling behavior of the CT with initial curvature, and it was shown that the initial curvature has a significant effect on the load in sinusoidal buckling and helical buckling [25,26].Subsequently, Qin and Gao investigated the buckling deformation of the tubular string with two initial configurations, sine and cosine, and found that the critical load values of sinusoidal buckling for the two initial configurations are basically the same [27].Zheng and Adnan assumed the buckling deformation of the tubular string when initially helical [28].Huang and Gao considered the effects of having both initial curvature and joints on the critical buckling load of the tubular column [29][30][31].
In this paper, a new method is used to intensively study the sinusoidal post-buckling deformation of the coiled tubing with residual bending in horizontal well.Firstly, a mathematical model of the bending deformation of a horizontal section of the coiled tubing, including a dynamic term, was established.The equations are reasonably simplified, and the separation variable is introduced to obtain an analytical solution for sinusoidal buckling with a time term.Secondly, when the coiled tubing with initial residual bending undergoes sinusoidal postbuckling deformation, the variation rule of the separation variable is obtained by using the energy method.The formulae for the compressive displacement and axial load of the coiled tubing containing the separation variable and time term are derived.Finally, contact force change rule analyzed for coiled tubing undergoing sinusoidal post-buckling with residual bending and motion speeds.

Mathematical model
Before constructing the mathematical model, the following assumptions must be made to simplify the analysis process and ensure the accuracy of the results.
1.The CT cannot be compressed; 2. The mechanical effects of external factors such as downhole fluids on the CT are ignored; 3. The cross sections of the wellbore and the CT are circular; 4. The CT is always in contact with the well wall when deformation occurs; 5.The sinusoidal buckling of the CT satisfies the small deformation theory; 6.The effect of the CT rotation and torsion is ignored.
7. Both ends of the CT are connected in a hinged mode.
The establishment of the coordinate system is shown in Fig 1 .whereθ is the angular displacement of the CT when buckling deformation occurs, and N is the contact force between the well wall and the CT.In the deformation process of the horizontal coiled tubing, any microsegment of the CT is taken as the research object, and the motion balance equation of the CT per unit length is established as follows: According to the momentum theorem: where V is the internal force vector of the microsegment, m is the mass of the microsegment, F is the vector of external forces per unit length of the CT, Δs is the microsegment arc length, v is the motion velocity vector of the microsegment center point O 2 .The moment of momentum theorem of the microsegment for wellbore center O 1 can be expressed as: where r o is the vector diameter from the origin O to the center of the wellbore O 1 , M is the external torque of the external force on the CT per unit length to the center O 2 of the CT, ψ is the internal moment on the microsegment of the CT, r is the vector diameter from the origin O to the center of the element O 2 , H is the moment of momentum per unit length of the CT concerning wellbore center O 1 .
From Eq (1)-( 3), and ignoring small quantities, we can obtain Eq (4) as follows: Assuming that e is a tangential unit vector, the two ends of the equal sign in Eq (4) are divided by Δs and taken as infinitesimal, and then: The angular momentum of the microsegment moving along the well wall is given by Eq (6).
where I c is the moment of inertia per unit length of the CT around its own axis, ω is the angular velocity vector of the rotation of the CT on its axis, O is the angular velocity vector of the centroid of the unit length of the tubular string rotating in the XY-plane about the center of the wellbore.The vector parameters are represented by the subscripts X, Y, and Z, corresponding to components along the X-axis, Y-axis, and Z-axis, respectively.When the CT is buckled and deformed, Gao and Mitchell believe that the movement direction of the contact point between the CT and the well wall can be divided into two directions: the movement direction along the Z-axis and the well wall rotation along the XY-plane [16,33].It is assumed that the dynamic friction coefficient along the Z-axis is μ z and the dynamic friction coefficient in the XY-plane is μ xy .The expressions of the friction vector f and the external force vector F per unit length of the CT string in continuous contact with the well wall are as follows respectively: where sign (θ) = 1 when the CT microsegment rotates clockwise, and sign(θ) = -1 when the CT microsegment rotates counterclockwise.
Assuming M is the external torque of the external force on the unit-length tubular string to the center O 2 of the CT, and r o2c is the vector from the contact point C to the center O 2 of the tubular string microsegment, as shown in Fig 1 .From Eq (5), ( 6) and ( 9), it can be obtained: If the radius is r, and r p is the CT radius, and r s is the wellbore radius, then: r = r s -r p .The following conditions are satisfied.
The symbol ()' is the derivative of the arc length s, and the component of the external force vector on the Z-axis is F. Eq (1)-( 6) are substituted into the constitutive equation of Gao for derivation [21], and ignoring torque, we can then substitute Eq ( 7)-( 11) and use the Lagrange formula to expand the derivation and obtain: 3y 000 y 03 r 3 cosy þ 9y 02 y @ 2 r 3 cosy À 3y 04 y @ r 3 siny À 2rz 0 z @ y 03 siny þ6rz 0 z @ y 0 y @ cosy þ 2rz 0 z @ y 000 siny À rz @ z 000 y 0 siny À rz 02 y 04 cosy þ3rz 02 y @ 2 cosy þ 4rz 02 y 0 y 000 cosy À 6rz 02 y 02 y @ siny þ rz 02 y @@ siny À rz 0 z 000 y 02 cosy À rz 0 z 000 y @ siny À rz 0 z @@ y 0 siny 3y 000 y 03 r 3 siny þ 9y 02 y @ 2 r 3 siny þ 3y 04 y @ r 3 cosy þ 2rz 0 z @ y 03 cosy þ6rz 0 z @ y 0 y @ siny À 2rz 0 z @ y 000 cosy þ rz @ z 000 y 0 cosy À rz 02 y 04 siny þ3rz 02 y @ 2 siny þ 4rz 02 y 0 y 000 siny þ 6rz 02 y 02 y @ cosy À rz 02 y @@ cosy À rz 0 z 000 y 02 siny þ rz 0 z 000 y @ cosy þ rz 0 z @@ y 0 cosy where: where E is the elastic modulus, and I is the moment of inertia.According to the hypothesis of Li [34], ds�dz.From Eqs ( 12) and ( 13), we can get the fourth-order nonlinear partial differential equations including the dynamic effect when the coiled tubing in the horizontal well is compressed and bent by axial load: Eq (14) are the same as the conclusion of Li [35].If the CT is stationary, equations become the form of inclination angle β = 0, according to Mitchell's conclusion.Because the CT is in continuous contact with the well wall, the contact force changes little with the axial coordinate z, and the friction coefficient μ z is less than 1.Additionally, because the contact force on the well wall at the beginning of the buckling deformation of the CT is approximately equal to the gravity of the CT, the values of μ z Nr p θ 0 and μ z r p N 0 in ( 14) are very small.It shows that the influence of axial friction on the tubular string deformation is far less than that of radial friction at the beginning of the tubular string deformation, which is the same as that of Gao and Miska [23,33].

Dimensionless
The following variables are defined: The dimensionless parameters are: where q is gravity per unit length of the CT.Because the initial deformation is small, for the convenience of solving, ignoring the influence of friction and high-order small quantity [22], Eqs ( 15) is a non-dimensional partial differential equation set for the deformation of horizontal coiled tubing, taking into account dynamic effects.

Energy analysis
The two ends of the horizontal coiled tubing are connected as hinged, then the boundary conditions are: In the boundary condition of the equation, the axial dimensionless length of the horizontal string after deformation is � L. Assuming that the starting end of the coiled tubing is hinged and the tail is moving, an analytical solution to the Eq (15) is obtained by substituting boundary conditions and introducing a separation variable, as shown in Eq (16).
where λ is the separation variable, C 1 and C 2 are arbitrary constants, and β = 2r4 / I.The necessary and sufficient conditions to satisfy the above equation are When λ = 1, Eq (17) becomes the Euler critical load hinged at both ends.When the CT is in the static state or uniform motion state, it satisfies: @ 2 θ / @ � t 2 = 0.The definition of the separation variable, λ = 0, satisfies this requirement.The expression for axial load excludes the velocity term, which is the same as the conclusion of Gao [22].According to the conclusion of Haberman [37], when solving the equation with a small angle, that is, sin (θ) � θ, the separation variable must satisfy λ < 0 when the CT is sinusoidal buckling.Eq (16) can be simplified as: where: For the CT with initial residual bending can be expressed by Eq (19).
where A 0 is the initial deflection of the CT, and m is the initial half-wave number, and L 0 is the initial axial length of the CT.The conditions are satisfied at this time.
The x 0 and y 0 are the horizontal and vertical coordinates values corresponding to the initial state of the CT, respectively, and θ 0 is the initial angular displacement of the cross-section.
The energy of the CT with buckling deformation are as follows [10,38,39].
where ΔU is the bending strain energy, ΔW q is the gravitational potential energy, ΔW F is the axial load work.Because the integral term appears to transcend functions sin(Asin θ) and cos (Asin θ).To simplify the calculation, the integral term is expanded using a power series as shown in Eq (23).
Using the potential energy principle, we get: By substituting Eq (20)-( 23) respectively, and dimensionless, we can get According to the principle of the minimum potential energy of a system, there is Due to the arbitrariness of δA n , the relationship between the dimensionless axial load and the amplitude A n of the horizontal coiled tubing is obtained after @P / @A n = 0 is satisfied Eq (25).
Taking the derivative of � L at both ends of the Eq (25) and setting it equal to zero, we can get By substituting Eq (26) into (25), the expression of minimum dimensionless axial load required for buckling deformation of a horizontal coiled tubing with initial deflection A 0 and initial half-wave number m can be obtained.
If the initial mode of the CT is a straight-line shape, then there is A 0 = 0. Eq (27) becomes: At this point, the CT just deformed can be considered A n = 0, the axial load becomes: This result is the same as the conclusion of Palsay, and the conclusion in Miska when α = 90˚ [3,32].
The dimensionless expression for the compression displacement of a horizontal coiled tubing is as follows: Substituting Eq (26) into (28) yields an expression for the amplitude A n and the dimensionless compression displacement � L z as shown in Eq (29).
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 9ðC � where:

Numerical analysis
Assuming the outsider diameter of the coiled tubing is 2 inches, the wall thickness of the CT is 0.165 inches, the radius of the wellbore is 3.5 inches, the density of the coiled tubing is 489 lb/ ft 3 , and the modulus of elasticity is E = 2.988×10 7 PSI.The relationship between the compression displacement and the generated amplitude at different initial amplitudes can be obtained from Eq (29), as shown in  the same, and the increase of A n under the same displacement of compression is the same.From Eq (29) we can see that the initial variable Г is the main factor influencing the amplitude A n .Changing the initial half-wave number and changing the initial amplitude both change the Г value.If the change of Г value is small, the effect on the change of amplitude A n is small.Fig 3 (A 0 = 0.5, m = 2) shows the amplitude change curves of columns with different initial lengths during compression.It can be seen that the longer the initial length of the CT, the smaller the amplitude under compression for the same displacement.This is because the longer CT generates more half-waves in compression [35], causing the amplitude to decrease.
By substituting Eq (29) into (25), we can obtain the dimensionless axial load � F versus compression displacement � L z required for bending of horizontal coiled tubing with different initial half-wave numbers, as shown in Fig 4(A) ( � L 0 = 100).The axial load necessary for the CT to experience buckling deformation increases with the number of half-waves, however, the increase is relatively minor.It is evident that the axial load is less affected by the change in the initial half-wave number for smaller Г values.If we continue to increase the value of Г and shorten the dimensionless length of the CT to 50, the relationship between the dimensionless axial load and the compression displacement is shown in Fig 4(B).
The initial half-wave number increases, the less axial load is required to compress the same displacement.However, the critical load for sinusoidal buckling deformation of the CT with initial deflection exceeds the critical load value for sinusoidal buckling of the initially straight CT.The relationship between axial load and compression displacement at varying initial amplitudes as the Г value increase is depicted in  initial half-wave number m = 40.For the CT with m = 50, the dimensionless compression displacement is approximately 0.11 when the axial load is extreme.This means that the larger the initial half-wave number is, the earlier the axial load extreme occurs.
Actually, Eq (27) takes the derivative of the amplitude A n at both ends and sets it equal to zero, we can get: Since the CT has residual bending, the solution for the amplitude is obtained as follows: It can be seen that the axial load � F will have an extreme point when Г is taken at the appropriate value as shown in Eq (30).
Li shows that the critical point where sinusoidal buckling transforms into sinusoidal-helical buckling is the extreme point of Eq (30) [35].As the value of Г increases, the axial load will appear earlier at the critical point, and the sinusoidal-helical shape will appear earlier as well.As compression continues, the angular displacement θ increases, and causing the error resulting from Eq (14) with approximate values to becomes larger.The curve behind the extreme point of the axial load in Fig 4(C) has a large error and is not practically relevant.The larger the length of the CT with initial deflection, the smaller the axial load is required to compress the same displacement when buckling.This conclusion is the same as the rule of sinusoidal buckling deformation of horizontal coiled tubing with straight initial shape considered by Li [35].

Influence of separation variable
Substituting Eqs ( 25), ( 26) and ( 29) into Eq (17) gives the relationship between the separation variable and the amplitude.L 0 = 10, the curve of compression displacement versus separation constant during compression of dimensionless displacement to 0.2 is shown in Fig 9 .According to the analytical results of Timoshenko [36], the deformation of the CT occurs due to the presence of the basic reaction force, so the separation variable is satisfied: λ � 1.If λ > 1, the contact force of the well wall on the buckling coiled tubing is negative, and at this time, the CT is separated from the contact with the well wall, and even more severe plastic deformation occurs, and it is obvious that the calculation using the formula is no longer meaningful.The larger the value of the initial parameter Г, the earlier the extremes of λ will occur.

Numerical calculation of contact force
The expression for the dimensionless contact force � N on the well wall when sinusoidal buckling occurs in coiled tubing with residual bending can be obtawined from the coiled tubing deformation Eq (16) and the coiled tubing initial configuration Eq (19) by substituting them into Eq (15) as shown in Eq (31).The symbol definitions in the expressions are consistent with the settings in the paper.
ffi ffi ffi ffi ffi ffi ffi  It can be seen that the contact force of the CT on the well wall changes periodically with the change in length and is greater than zero.This result is the same as the conclusion of Gao [21].By comparing the contact force curve of different lengths of the tubular string on the well wall, it can be seen that other parameters being equal, the larger the length of the CT with residual bending, the larger the maximum contact force on the wellbore under the same displacement, while the minimum contact force is almost not affected by the length of the CT.The maximum value of the contact force on the well wall of a buckling tubular column increases as the initial amplitude increases, and the increase is greater.When A 0 = 0.4, � N = 4.57, when A 0 = 0.6, � N = 22.54, when A 0 = 0.8, � N = 4.57.When A 0 = 0.8, the maximum contact force is reached � N max = 55.9.Therefore, the maximum contact force on the well wall when the CT with residual bending is buckled is greater than that when the CT with a straight initial The greater the compression displacement, the smaller the value of the maximum contact force of the CT on the well wall and the larger the value of the minimum contact force.The compression displacement becomes larger, the axial load required necessarily increases, and the contact force on the CT located at the bottom of the wellbore increases.However, the CT tends to return to its initial form at this point, causing the contact force of the sinusoidal peak to gradually decrease.The contact force period decreases over time, indicating a reduction in the number of sinusoidal half-waves generated by the buckled coiled tubing as compression progresses.Eventually, the CT transforms into a mixed sinusoidal-helical buckling deformation pattern.velocity � v = 0.08, the maximum value of dimensionless contact force increases to � N = 2619.If the velocity is further increased to � v = 0.1, the maximum value of dimensionless contact force reaches � N = 4092.Therefore, if the speed at which the CT is lowered exceeds a certain value during actual production, the contact force of the buckling coiled tubing on the well wall increases rapidly.This is one of the reasons why the coiled tubing lowering speed cannot be too fast.

Conclusions
This paper presents a generalized mathematical model for the buckling deformation of a horizontal coiled tubing during its motion.The model is established using the momentum theorem and the momentum moment theorem.The mathematical model includes the time term for the sinusoidal buckling deformation of a coiled tubing with residual bending.The equations are simplified and a separation variable is introduced.The residual bending of the coiled tubing is considered a sinusoidal shape.The amplitude, axial load, and contact force variation of the coiled tubing after sinusoidal buckling are analyzed using the energy method.This lays the foundation for further research on the process of the horizontal coiled tubing with initial deflection transforming from sinusoidal buckling to helical buckling, and the following conclusions are obtained.
1. Theoretical analysis was conducted on the emergence conditions of sinusoidal buckling of horizontal coiled tubing with initial deflection to mixed sinusoidal-helical buckling critical point.The results indicate that as the initial parameter Γ of the coiled tubing increases, the transition from sinusoidal buckling to sinusoidal-helical buckling critical point occurs earlier during compression.If the initial parameter Γ is very small, the sinusoidal-helical buckling deformation of the coiled tubing is not easy to occur.
2. The effects of compression displacement on the amplitude, axial load, and separation variable of the coiled tubing with different initial amplitudes and initial half-wave numbers are analyzed.The deformation press is similar to that calculated by Wu when the Γ value is small, and the initial half-wave number and initial amplitude have little influence on the deformation of the coiled tubing.As the value of Γ increase, the critical load of sinusoidal buckling of the coiled tubing also increases with the initial amplitude.However, the change of the initial half-wave number has little effect on the critical load of sinusoidal buckling.
3. A mathematical model was developed to calculate the contact force of a horizontal coiled tubing with residual bending subjected to sinusoidal post-buckling, which includes a time term.The contact force changes on the well wall during sinusoidal post-buckling of the coiled tubing under different conditions were calculated.According to the calculation results, when the coiled tubing with residual bending undergoes buckling deformation, the maximum contact force on the well wall is greater than that calculated by Gao.The longer the coiled tubing, the greater the maximum contact force on the well wall at the same distance of compression.The more the initial half-wave number, the greater the maximum contact force.The larger the initial amplitude, the greater the maximum value of contact force.However, the minimum value of contact force under the above conditions is not affected.Finally, it was proven theoretically that the lowering speed of the coiled tubing has a significant impact on the contact force.and the faster the lowering speed, the greater the contact force.
Mathematically, the set of differential equations for the buckling deformation of the coiled tubing is a system of high-order strongly nonlinear coupled system of partial differential equations, so the process of solving the equations uses a more idealized case, ignoring some boundary conditions such as friction, torque, and some small values.Although the compression distance at which the CT undergoes sinusoidal buckling is small, its deformation increases when the CT continues to be compressed, which may result in large errors if the small angle approximation in the paper continues to be used for analytical calculations.Therefore, the analysis process in this paper is limited to the process of sinusoidal buckling deformation of the CT.
The future research will focus on the deformation state of the CT when the spiral buckling occurs using the mathematical model established in this paper, the effect of increasing the friction, torque, and other boundary conditions on the mathematical model, and how the coiled tubing is converted from sinusoidal buckling to spiral buckling in the process of motion.

Fig 2 (
A 0 = 0.5, � L 0 = 100).The change of the initial half-wave number has little effect on the change of amplitude A n , which is basically consistent with the results of Wu[40].The slopes of the amplitudes are all

Fig 5 (
A) illustrates the relationship between axial load and compression displacement for sinusoidal buckling of the CT with various initial amplitudes, where the initial half-wave number m = 2 and dimensionless length � L 0 = 100.Fig 6 shows the same initial conditions (m = 2, A 0 = 0.2), the relationship between axial load and compression displacement of different initial lengths of the CT during compression.
Fig 7(A)-7(C) show the relationship between separation variable

Fig 11 (
A) and 11(B) show the relationship curves of the contact force of the CT on the well wall when compressed with the same displacement under different initial amplitudes A 0 ( � L 0 = 50, m = 5, � L z = 0.005, � t = 1) and different initial half-wave numbers m ( � L 0 = 50, A 0 = 2, � L z = 0.005, � t = 1) at the same initial length, respectively.

Fig 10 .
Fig 10.Contact force curves of the CT with different lengths on the horizontal well wall.https://doi.org/10.1371/journal.pone.0301610.g010 Fig 13 shows the variation of the contact force at a compression of 0.01 dimensionless displacement for the CT with dimensionless length � L 0 = 200, initial amplitude A 0 = 0.2, and initial half-wave number m = 2.The dimensionless velocities � v at the end is taken as 0.02, 0.04, 0.06, and 0.08, respectively.Fig 13 shows that the greater the coiled tubing lowering speed, the greater the maximum contact force generated by the deformation of the coiled tubing on the well wall.The minimum value remains largely unchanged and is estimated as � N = 0.8171.The maximum value of dimensionless contact force is � N = 165.3when � v = 0.02.However, when the dimensionless